Library
library(AICcmodavg)
library(ape)
library(caper)
library(geiger)
library(kableExtra)
library(MuMIn)
library(nlme)
library(phytools)
library(plotrix)
library(shape)
library(car)
library(magick)Session information
R.version
_
platform x86_64-apple-darwin17.0
arch x86_64
os darwin17.0
system x86_64, darwin17.0
status
major 4
minor 2.2
year 2022
month 10
day 31
svn rev 83211
language R
version.string R version 4.2.2 (2022-10-31)
nickname Innocent and Trusting
sessionInfo()
R version 4.2.2 (2022-10-31)
Platform: x86_64-apple-darwin17.0 (64-bit)
Running under: macOS Catalina 10.15.7
Matrix products: default
BLAS: /Library/Frameworks/R.framework/Versions/4.2/Resources/lib/libRblas.0.dylib
LAPACK: /Library/Frameworks/R.framework/Versions/4.2/Resources/lib/libRlapack.dylib
locale:
[1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
attached base packages:
[1] stats graphics grDevices utils datasets methods base
other attached packages:
[1] magick_2.7.4 car_3.1-1 carData_3.0-5 shape_1.4.6
[5] plotrix_3.8-2 phytools_1.5-12 maps_3.4.1 nlme_3.1-162
[9] MuMIn_1.47.5 kableExtra_1.3.4 geiger_2.0.10 caper_1.0.1
[13] mvtnorm_1.1-3 MASS_7.3-58.3 ape_5.7-1 AICcmodavg_2.3-2
[17] rmarkdown_2.20.1
loaded via a namespace (and not attached):
[1] Rcpp_1.0.10 subplex_1.8 svglite_2.1.1
[4] lattice_0.20-45 digest_0.6.31 unmarked_1.2.5
[7] foreach_1.5.2 R6_2.5.1 plyr_1.8.8
[10] stats4_4.2.2 evaluate_0.20 coda_0.19-4
[13] httr_1.4.5 rlang_1.1.0 rstudioapi_0.14
[16] jquerylib_0.1.4 phangorn_2.11.1 Matrix_1.5-3
[19] combinat_0.0-8 splines_4.2.2 webshot_0.5.4
[22] stringr_1.5.0 munsell_0.5.0 igraph_1.4.1
[25] compiler_4.2.2 numDeriv_2016.8-1.1 xfun_0.37
[28] systemfonts_1.0.4 pkgconfig_2.0.3 mnormt_2.1.1
[31] optimParallel_1.0-2 htmltools_0.5.4 expm_0.999-7
[34] quadprog_1.5-8 codetools_0.2-19 viridisLite_0.4.1
[37] grid_4.2.2 lifecycle_1.0.3 jsonlite_1.8.4
[40] xtable_1.8-4 magrittr_2.0.3 scales_1.2.1
[43] stringi_1.7.12 cli_3.6.0 cachem_1.0.7
[46] pbapply_1.7-0 scatterplot3d_0.3-43 doParallel_1.0.17
[49] xml2_1.3.3 bslib_0.4.2 generics_0.1.3
[52] fastmatch_1.1-3 deSolve_1.35 iterators_1.0.14
[55] tools_4.2.2 glue_1.6.2 abind_1.4-5
[58] parallel_4.2.2 fastmap_1.1.1 survival_3.5-5
[61] yaml_2.3.7 colorspace_2.1-0 rvest_1.0.3
[64] VGAM_1.1-8 clusterGeneration_1.3.7 knitr_1.42
[67] sass_0.4.5Consider only sit-and-wait and widely-foraging lizards
Almost a year ago, my mentors and I submitted a paper to Proceedings of the Royal Society B. In this paper, we investigated the relationship between foraging behavior and reproducitve effort in lizards (see here). One of the most important comments I got from a reviwer was related to the scatter in the data. The reviwer stated that the variance of the data was very large, preventing us from making robust conclusions. In response to this concern, I would like to show here one way to deal with this issue:
As observed in figure 3a-b of the paper, the majority of the species are clustered at small sizes along the x-axis, while a few others are grouped at large sizes and seem to be “extreme cases”. Some of the latter cases are “mixed-foraging” species that can behave either as a “widely-foraging” or “sit-and-wait” species. If I remove the “mixed” category and make the species behave randomly as “widely foragers” or “sit-and-wait forager”, I could recalculate the significance of the interaction to determine whether the results hold or not. If I do so multiple time (>1000), I can calculate the proportion of time that “widely-foraging” species have greater reproductive output than “sit-and-wait” species, as indicated by the original results. Let’s do it:
## Considering only sit-and-wait and widely-foraging lizards
set.seed(94)
options(scipen = 999)
active <- dat_full[dat_full$foraging.mode == "widely foraging", ]
ambush <- dat_full[dat_full$foraging.mode == "sit and wait", ]
mixed <- dat_full[dat_full$foraging.mode == "mixed", ]
savep <- rep(NA, 1000)
savei <- rep(NA, 1000)
for(i in 1:1000){
randomizing <- sample(rep(c("widely foraging", "sit and wait"), each = length(mixed$foraging.mode)), size = length(mixed$foraging.mode))
mixed$foraging.mode <- randomizing
merg <- rbind(active, ambush, mixed)
## model
mod <- gls(rep_out ~ M_females*foraging.mode, correlation = corPagel(value = 0, phy = ctree_full, form = ~Binomial), data = merg, method = "ML")
## Confidence intervals
## sit and wait
sw <- merg[merg$foraging.mode == "sit and wait", ]
check_sw <- name.check(tree, sw)
rm_phy_sw <- check_sw$tree_not_data
rm_dat_sw <- check_sw$data_not_tree
ctree_sw <- drop.tip(tree, rm_phy_sw)
SSX.sw <- sum(round((sw$M_females - mean(sw$M_females))^2), 2)
s2.sw <- var(sw$rep_out)
n.sw <- length(sw$rep_out)
x.sw <- seq(min(sw$M_females), max(sw$M_females), 0.5)
m.x.sw <- mean(round(sw$M_females, 1))
se.sw <- sqrt(s2.sw*((1/n.sw) + (((x.sw - m.x.sw)^2)/SSX.sw)))
is.sw <- qt(0.975, df = n.sw - 2)
ii.sw <- qt(0.025, df = n.sw - 2)
ic.s.sw <- se.sw*is.sw
ic.i.sw <- se.sw*ii.sw
upper.i.sw <- ((coef(mod)[1] + coef(mod)[3]) + (coef(mod)[2] + coef(mod)[4])*x.sw) + ic.s.sw
lower.i.sw <- ((coef(mod)[1] + coef(mod)[3]) + (coef(mod)[2] + coef(mod)[4])*x.sw) + ic.i.sw
## widely foraging
wf <- merg[merg$foraging.mode == "widely foraging", ]
check_wf <- name.check(tree, wf)
rm_phy_wf <- check_wf$tree_not_data
rm_dat_wf <- check_wf$data_not_tree
ctree_wf <- drop.tip(tree, rm_phy_wf)
SSX.wf <- sum(round((wf$M_females - mean(wf$M_females))^2), 2)
s2.wf <- var(wf$rep_out)
n.wf <- length(wf$rep_out)
x.wf <- seq(min(wf$M_females), max(wf$M_females), 0.5)
m.x.wf <- mean(round(wf$M_females,1))
se.wf <- sqrt(s2.wf*((1/n.wf) + (((x.wf - m.x.wf)^2)/SSX.wf)))
is.wf <- qt(0.975, df = n.wf - 2)
ii.wf <- qt(0.025, df = n.wf - 2)
ic.s.wf <- se.wf*is.wf
ic.i.wf <- se.wf*ii.wf
upper.i <- (coef(mod)[1] + coef(mod)[2]*x.wf) + ic.s.wf
lower.i <- (coef(mod)[1] + coef(mod)[2]*x.wf) + ic.i.wf
png(paste("figure", i, ".png", sep = ""), width = 7, height = 7, units = "in", res = 300)
plot(merg$rep_out ~ merg$M_females, type = "p", xlab = " ", ylab = "Hatchling mass x clutch size (g)", cex.lab = 1.2, pch = 21, bg = c("purple", "skyblue")[as.numeric(merg$foraging.mode)], col = c("purple", "skyblue")[as.numeric(merg$foraging.mode)], las = 1)
polygon(c(rev(x.sw), x.sw), c(rev(lower.i.sw), upper.i.sw), border = FALSE, col = rgb(0, 0, 1, alpha = 0.2))
lines(x = x.sw, y = ((coef(mod)[1] + coef(mod)[3]) + (coef(mod)[2] + coef(mod)[4])*x.sw), col = "skyblue", lwd = 2)
polygon(c(rev(x.wf), x.wf), c(rev(lower.i), upper.i), border = FALSE, col = rgb(0.8, 0, 0.7, alpha = 0.3))
lines(x = x.wf, y = coef(mod)[1] + coef(mod)[2]*x.wf, col = "purple", lwd = 2)
text(x = 7.7, y = 30, substitute(paste(italic("Interaction \n "))), family = "serif")
text(x = 7.7, y = 30, round(coef(mod)[4], 3), family = "serif")
text(x = 7.7, y = 27, paste("p = ", round(anova(mod)[4, 3], 3), sep = " "), family = "serif")
mtext(expression(hat(M) (g)), side = 1, line = 4, cex = 1)
legend("topleft", legend = c("sit-and-wait foraging", "widely foraging"), lty = 1, bty = "n", pch = c(16, 16), bg = c("black", "black"), col = c("skyblue", "purple"), cex = 1)
dev.off()
savep[i] <- round(anova(mod)[4, 3], 3)
savei[i] <- round(coef(mod)[4], 3)
}
## Estimate proportion of time that the p-value of the interaction was significant
mean(savep < 0.05)
[1] "0.574"
## Mean affect size of the interaction
mean(savei)
[1] "-0.259"
## list file names and read in
imgs <- list.files(path = "/Users/dpadil10/Dropbox (ASU)/gif/gif/", full.names = TRUE)
img_list <- lapply(imgs, image_read)
## join the images together
img_joined <- image_join(img_list)
## animate at 2 frames per second
img_animated <- image_animate(img_joined, fps = 2)
img_animated